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Combustion  Research  Facility 
Livermore,  California  94551-0960 


Abstract 

This  paper  provides  a perspective  on  Large  Eddy  Simulation  (LES)  and  its  application  to  liquid  rocket  injection, 
mixing  and  combustion  processes.  Simulating  these  processes  involves  a variety  of  challenges  which  include  all  of 
the  classical  closure  problems  inherent  to  the  treatment  of  turbulence  and  combustion,  and  a unique  set  of  problems 
imposed  by  the  introduction  of  thermodynamic  nonidealities  and  transport  anomalies.  Emphasis  is  placed  on  1)  the 
fundamental  issues  and  limiting  extremes,  2)  the  theoretical  and  numerical  framework  developed  to  handle  these 
difficulties,  and  3)  a series  of  results  which  give  insight  into  the  intricate  nature  of  the  problem  and  current  state  of 
the  art  with  respect  to  LES.  The  discussion  is  framed  in  the  context  of  the  three  workshop  test  cases,  with  conclusions 
drawn  accordingly. 


Introduction 

Simulating  injection,  mixing  and  combustion  processes  in  cryogenic  rocket  engines  poses  a variety  of  challenges  which 
include  all  of  the  classical  closure  problems  inherent  to  the  treatment  of  combustion,  and  a unique  set  of  problems 
imposed  by  the  introduction  of  thermodynamic  nonidealities  and  transport  anomalies.  Flow  conditions  within  the 
chamber  are  inherently  turbulent,  and  combined  demands  associated  with  performance  and  heat  transfer  typically 
result  in  the  specification  of  operating  pressures  and  temperatures  that  produce  local  transcritical 1 and  supercritical 
conditions. 

From  the  classical  point  of  view,  reacting  multiphase  flows  introduce  the  complicating  factors  of  chemical  kinetics, 
highly  nonlinear  source  terms,  and  a variety  of  subgrid-scale  (sgs)  velocity  and  scalar  mixing  interactions.  Flow  field 
evolution  is  affected  by  compressibility  effects  (volumetric  changes  induced  by  changes  in  pressure)  and  variable 
inertia  effects  (volumetric  changes  induced  by  variable  composition  and/or  heat  addition).  The  resultant  coupling 
yields  an  array  of  fluid  dynamic,  thermochemical,  thermodynamic  and  transport  processes  which  are  dominated  by 
widely  disparate  time  and  length  scales. 

The  situation  becomes  more  complex  at  elevated  pressures  due  to  the  inherent  decrease  in  turbulence  scales  and 
difficulties  which  arise  as  fluid  states  approach  and/or  exceed  local  critical  conditions.  Near  the  critical  point,  propel- 
lant mixture  properties  begin  to  exhibit  liquid-like  densities,  gas-like  diffusivities  and  pressure-dependent  solubilities. 
Surface  tension  and  heat  of  vaporization  approach  zero,  and  the  isothermal  compressibility  and  constant  pressure  spe- 
cific heat  increase  significantly.  These  phenomena,  coupled  with  extreme  local  property  variations,  have  a significant 
impact  on  the  evolutionary  dynamics  of  a given  system. 

This  paper  presents  a perspective  on  the  Large  Eddy  Simulation  (LES)  technique  and  its  application  to  the  injection, 
mixing  and  combustion  processes  described  above.  After  establishing  the  key  phenomenological  trends  and  flow 
characteristics,  the  implications,  modeling  options,  and  tradeoffs  are  outlined  and  the  general  requirements  for  LES 
are  discussed  and  contrasted  with  the  Reynolds-Averaged  Navier-Stokes  (RANS)  approach.  A series  of  case  studies 
are  then  presented  which  highlight  various  aspects  of  the  problem  and  the  utility  of  LES  as  a fundamental  tool.  The 
discussion  is  framed  in  the  context  of  the  three  workshop  test  cases,  with  conclusions  drawn  accordingly. 

1 A liquid  propellant  at  subcritical  temperature  in  a high-pressure  supercritical  environment. 
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Figure  I;  Reacting  shear-coaxial  liquid-oxygen-hydrogen  injector  operating  at  1.5  MPa  (15  atm).  From  Mayer  and 
Tamura  (1996).  Used  with  permission. 


Figure  2:  Reacting  shear-coaxial  liquid-oxygen-hydrogen  injector  operating  at  4.5  MPa  (44  atm).  From  Mayer  and 
Tamura  (1996).  Used  with  permission. 


Phenomenological  Trends 

Recent  experiments  have  provided  a much  clearer  understanding  of  the  phenomenological  conditions  which  exist  as  a 
function  of  chamber  pressure.  Depending  on  the  injector  type,  fluid  properties,  and  flow  characteristics,  two  limiting 
extremes  may  exist.  At  subcritical  chamber  pressures,  injected  liquid  jets  undergo  the  classical  cascade  of  processes 
associated  with  atomization.  Dynamic  forces  and  surface  tension  promote  the  formation  of  a heterogeneous  spray 
which  evolves  continuously.  Spray  flames  form  as  a consequence  which  are  lifted  away  from  the  injector  face  in 
a manner  consistent  with  the  combustion  mechanisms  exhibited  by  local  drop  clusters.  When  chamber  pressures 
approach  or  exceed  the  critical  pressure  of  a particular  propellant,  however,  injected  liquid  jets  undergo  a transcriticai 
change  of  state  as  interfacial  fluid  temperatures  rise  above  the  critical  temperature  of  the  local  mixture.  For  this 
situation,  diminished  intermolecular  forces  promote  diffusion  dominated  processes  prior  to  atomization  and  respective 
jets  vaporize  forming  a continuous  fluid  in  the  presence  of  exceedingly  large  gradients.  Well  mixed  diffusion  flames 
evolve  as  a consequence  which  are  anchored  by  small  but  intensive  recirculation  zones  generated  by  the  sheariayers 
imposed  by  adjacent  propellant  streams. 

The  flow  visualization  studies  conducted  by  Mayer  and  Tamura  (1996)  illustrate  these  trends  for  the  case  of  a 
liquid-oxygen-gaseous-hydrogen  shear-coaxial  injector  element.  The  two  extremes  are  shown  in  Figs.  1 and  2, 
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Figure  3.  Near  injector  region  of  a reacting  liquid-oxygen— gaseous-hydrogen  shear-coaxial  injector,  flame  (top)  and 
corresponding  flow  field  (bottom).  Oxygen  and  hydrogen  velocities  are  30  and  300  m/s,  respectively,  oxygen  and 
hydrogen  injection  temperatures  are  100  K and  300  K,  oxygen  jet  diameter  is  1 mm,  chamber  pressure  is  4.5  MPa  (44 
atm).  From  Mayer  and  Tamura  (1996).  Used  with  permission. 


respectively.  Note  that  the  critical  pressure  and  temperature  of  oxygen  are  pc  = 5.04  MPa  (49.7  atm)  and  Tc  = 155 
K,  respectively.  The  critical  pressure  and  temperature  of  hydrogen  are  p c = 1.30  MPa  (12.8  atm),  Tc  = 33.2  K. 
When  liquid-oxygen  is  injected  at  low-subcritical  pressures  (Fig.  1)  atomization  occurs  forming  a distinct  spray  as 
described  above.  Ligaments  are  detached  from  the  jet  surface  forming  spherical  drops  which  subsequently  breakup 
and  vaporize.  As  the  chamber  pressure  approaches  the  thermodynamic  critical  pressure  of  the  liquid-oxygen  (Fig. 
2),  the  number  of  drops  present  diminishes.  Here,  the  injected  jet  exhibits  a pure  diffusion  mechanism  at  a pressure 
of  4.5  MPa , which  is  slightly  below  the  thermodynamic  critical  pressure  of  oxygen,  and  significantly  above  that  of 
hydrogen.  Experimental  results  have  revealed  that  flame  attachment  occurs  instantaneously  after  ignition  in  the  small 
but  intensive  recirculation  zone  which  forms  just  downstream  of  the  annular  post.  A well  mixed  diffusion  flame 
forms  within  this  region  producing  a wake  that  separates  the  oxygen  stream  from  the  hydrogen-rich  outer  flow.  The 
conditions  imposed  for  Test  Cases  RCM-2  and  RCM-3  are  phenomenologically  analogous  to  those  of  Figs.  1 and  2. 

Flow  Characteristics 

Simulating  either  of  the  two  extremes  described  above  with  either  LES  or  RANS  based  methods  requires  a detailed 
representation  of  the  broadband  turbulence  coupled  with  appropriate  multiphase,  thermochemical,  thermodynamic 
and  transport  models.  Modeling  subcritical  atomization  and  dense  spray  processes  similar  to  those  depicted  in  Fig.  1 
are  still  one  of  the  most  difficult  and  evasive  topics  of  research.  At  this  point  only  relatively  crude  highly  empirical 
models  exist.  Dilute  spray  models,  however,  are  more  prevalent  and  can  be  very  accurate  at  low  subcritical  pressures. 
Results  shown  in  subsequent  sections  will  demonstrate  the  effectiveness  of  LES  in  treating  dilute  spray  dynamics  in 
a configuration  which  eliminates  the  ambiguities  associated  with  atomization.  Simulating  the  transcritical  jet,  on  the 
other  hand,  does  not  require  use  of  an  atomization  model,  but  does  require  a detailed  representation  of  the  physical 
properties. 
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Because  the  fluid  is  much  denser,  the  broadband  turbulence  characteristics  which  must  be  considered  are  clearly 
evident  in  Fig.  2.  Figure  3 is  a visualization  which  illustrates  the  near  injector  region  of  this  case  in  the  vicinity  of 
the  liquid-oxygen  post.  The  mean  flame  characteristics  are  shown  on  the  top  of  the  figure,  and  the  corresponding  flow 
field  is  shown  on  the  bottom.  The  oxygen  and  hydrogen  velocities  for  this  case  are  30  and  300  m/s,  respectively,  the 
oxygen  and  hydrogen  injection  temperatures  are  100  K and  300  K , the  oxygen  jet  diameter  is  1 mm , and  the  mean 
chamber  pressure  is  4.5  MPa  (44  atm).  Figures  4 and  5 show  the  corresponding  thermophysical  behavior  of  of  oxygen 
and  hydrogen  over  the  regimes  of  interest.  Plots  of  density,  specific  heat,  viscosity,  and  thermal  conductivity  are  given 
on  the  interval  40  < T < 1000  K for  pressures  of  1, 10,  50,  100,  200,  and  400  atmospheres.  Note  that  at  1000  K and 
above,  both  oxygen  and  hydrogen  exhibit  ideal  gas  behavior  and  the  pressure  effect  is  negligible.  As  the  temperature 
is  decreased  below  1000  AT,  however,  significant  nonidealities  are  introduced,  with  property  variations  associated  with 
oxygen  producing  the  most  significant  effects. 

Figure  6 shows  the  trends  associated  with  the  kinematic  viscosity.  The  effect  of  pressure  on  this  quantity  is 
particularly  significant  and  has  a direct  impact  on  the  characteristic  scales  associated  with  the  turbulence  field.  For 
both  oxygen  and  hydrogen,  an  increase  in  pressure  from  1 to  100  atmospheres  results  in  a corresponding  reduction  in 
the  kinematic  viscosity  of  up  to  three  orders  of  magnitude.  This  implies  a three  order  of  magnitude  increase  in  the 
characteristic  Reynolds  number.  Based  on  Kolmogorov’s  universal  equilibrium  theory  (Tennekes  and  Lumley  1972, 
Hinze  1975),  the  order  of  magnitude  of  the  Kolmogorov  microscale,  denoted  here  as  rjt,  and  the  Taylor  microscale, 
denoted  here  as  A^,  are  related  to  the  Reynolds  number  by 

0) 

h n 

Here  the  Reynolds  number  is  defined  as  Ret  = qtltfv  where  qt  = y/2kt/3.  The  term  qt  represents  the  turbulence 
intensity,  kt  the  sgs  kinetic  energy,  and  It  the  integral  length  scale.  The  relations  given  by  Eq.  ( 1 ) indicate  that  a three 
order  of  magnitude  decrease  in  the  kinematic  viscosity  results  in  2,25  and  1.5  order  of  magnitude  decreases  in  the 
Kolmogorov  and  Taylor  microscales,  respectively.  These  reductions  have  a direct  impact  on  the  overall  grid  density 
required  to  resolve  key  processes. 

Figure  7 shows  the  trends  associated  with  the  effective  mass  diffusivity.  When  the  pressure  is  increased  from  1 
to  100  atmospheres,  both  oxygen  and  hydrogen  exhibit  a two  order  of  magnitude  decrease  in  the  mass  diffusion  rate 
over  the  full  range  of  temperatures  plotted.  Oxygen  exhibits  a decrease  of  up  to  4 orders  of  magnitude  at  temperatures 
below  the  critical  mixing  temperature.  The  diminished  mass  diffusion  rates  coupled  with  the  liquid-like  densities 
which  dominate  at  high  pressures  significantly  alter  the  coupling  associated  with  local  combustion  characteristics  in 
the  vicinity  of  the  liquid-oxygen  jet. 

Qualitative  analysis  of  Fig.  3 correlated  with  the  trends  shown  in  Figs.  4-7  suggests  that  there  are  at  least  seven 
fundamentally  important  flow  characteristics  which  must  be  accounted  for:  1 ) dense  near-critical  and  supercritical  fluid 
mixture  properties,  2)  transient  broadband  turbulent  mixing  over  a wide  range  of  scales,  3)  high  pressure  chemical 
kinetics,  4)  strong  multicomponent  property  gradients,  5)  dominant  preferential  diffusion  processes,  6)  anomalous 
multiphase  interfaces,  and  7)  geometrically  dominated  (wall-bounded)  three-dimensional  evolution.  Treating  this  set 
of  characteristics  represents  a minimal  requirement  for  any  simulation-based  or  modeled  treatment  of  the  flow.  In 
general,  the  intricate  multiple-time,  multiple-length  scale  coupling  must  be  resolved  (or  modeled)  to  represent  the 
physics.  A time-accurate  treatment  of  turbulence,  thermochemistry,  thermodynamics  and  transport  properties  are 
essential,  and  disparate  turbulence  and  molecular  transport  processes  must  be  treated  simultaneously.  Algorithm 
design  and  high-performance  massively  parallel  computing  are  also  essential  elements. 


Modeling  Options  and  Tradeoffs 

There  are  currently  three  basic  choices  with  regard  to  the  simulation  approach.  The  widely  used  Reynolds-Averaged 
Navier-Stokes  (RANS)  approximation.  Large  Eddy  Simulation  (LES)  and  Direct  Numerical  Simulation  (DNS).  The 
least  numerically  intensive  is  RANS.  For  this  approach,  all  turbulent  motions  are  modeled.  The  closure  is  empirical 
and  based  on  scaling  arguments  which  apply  only  in  the  time-averaged  limit.  In  general,  predictions  are  highly 
sensitive  to  models  and  model  constants,  and  respective  constants  must  be  adjusted  and  tuned  for  every  flow.  LES  is 
a much  more  numerically  intensive  methodology,  but  offers  a higher  degree  of  accuracy  in  return.  For  this  approach. 
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Figure  4:  Density  (compared  with  experimental  data  points  obtained  by  Vargaftik  (1975)),  specific  heat,  viscosity,  and 
thermal  conductivity  versus  temperature  over  the  interval  40  < T < 1000  and  pressures  of  1,  10,  50,  100,  200,  and 
400  atmospheres  for  pure  oxygen. 


Figure  5:  Density,  specific  heat,  viscosity,  and  thermal  conductivity  versus  temperature  over  the  interval  40  < T < 
1000  and  pressures  of  1, 10,  50, 100,  200,  and  400  atmospheres  for  pure  hydrogen. 
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Figure  6:  Kinematic  viscosity  of  pure  oxygen  (left)  and  pure  hydrogen  (right)  over  the  temperature  interval  40  < T < 
1000  and  pressures  of  1, 10,  50, 100, 200,  and  400  atmospheres. 


Figure  7:  Effective  mass  diffusivity  of  pure  oxygen  (left)  and  pure  hydrogen  (right)  over  the  temperature  interval 
40  < T < 1000  and  pressures  of  1, 10,  50, 100,  200,  and  400  atmospheres. 
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the  large  energetic  scales  are  resolved,  and  the  subgrid-scales  are  modeled.  In  contrast  to  RANS,  LES  closures  are 
time-accurate,  the  models  tend  to  be  more  universal,  and  it  is  not  necessary  to  adjust  constants  for  every  flow.  With 
the  appropriate  grid  constraints  in  place,  the  use  of  dynamic  modeling  eliminates  the  need  for  any  model  constants. 
The  enhanced  accuracy,  however,  comes  with  a much  stricter  set  of  algorithmic  requirements.  Similarly,  DNS  is  the 
most  numerically  intensive.  For  this  approach,  all  scales  are  resolved,  and  no  modeling  is  required,  but  the  method  is 
severely  CPU  limited. 

There  are  three  baseline  considerations  which  currently  dictate  the  selection  of  RANS  versus  LES  versus  DNS  as  a 
solution  method:  1)  the  time  required  to  get  a solution,  2)  the  accuracy  of  a solution,  and  3)  the  feasibility  of  obtaining 
a solution.  Performing  a DNS  which  meets  the  seven  criteria  listed  in  the  previous  section  is  clearly  not  feasible  at  this 
point  in  time.  Thus,  tradeoffs  typically  revolve  around  RANS  versus  LES,  with  DNS  being  used  as  a more  fundamental 
tool  for  studying  extremely  small-scale  phenomena  in  highly  idealized  domains.  This  in  itself  is  limiting,  and  care 
must  be  taken  to  insure  DNS  simulations  of  this  type  are  truly  relevant  to  the  flow  phenomena  of  interest.  LES,  by 
definition,  is  an  inherently  three-dimensional  simulation  methodology  (as  is  DNS),  and  LES  grid  requirements  are 
much  stricter  than  those  for  RANS.  Well  proportioned  LES  grids  are  typically  sized  a factor  of  four  coarser  in  each 
coordinate  direction  than  an  equivalent  well  sized  DNS  grid  with  the  equivalent  order  of  accuracy.  Examples  of  well 
sized  LES  grids  are  given  in  subsequent  sections.  RANS  is  by  far  the  fastest  solution  technique,  but  also  the  least 
accurate  since  the  entire  system  is  essentially  a model.  However,  the  speed  and  minimal  resources  with  which  RANS 
solutions  can  be  obtained  is  often  an  invaluable  and  necessary  engineering  tool.  A secondary  consideration  is  the  fact 
that  the  LES  methodology  imposes  a much  stricter  set  of  numerical  requirements  and  constraints. 


General  Requirements  for  LES 

Improvements  in  computational  speed  and  capacity  over  the  past  several  years  has  made  the  application  of  LES  fea- 
sible for  increasingly  complex  flows.  This  method  has  now  been  used  successfully  as  both  as  a complementary  tool 
for  understanding  turbulence  and  for  modeling  the  effects  of  turbulence  in  a variety  of  engineering  applications.  With 
the  advent  of  massively  parallel  computer  hardware,  LES  now  provides  a means  to  study  coupled  combustion,  trans- 
port and  multiphase  processes  in  parameter  spaces  that  are  unattainable  using  direct  numerical  simulation  (DNS) 
techniques,  with  a degree  of  fidelity  that  can  be  far  more  accurate  than  other  conventional  methods. 

Modeling  Issues 

During  the  past  decade,  considerable  progress  has  been  made  in  LES.  Early  works  relied  heavily  on  the  use  of  the 
Smagorinsky  eddy  viscosity  model  (Smagorinsky  1963).  A key  breakthrough  in  the  field  of  sgs  modeling  resulted 
from  the  introduction  of  the  dynamic  modeling  procedure  (Germano,  Piomelli,  Moin  and  Cabot  1991,  Moin,  Squires, 
Cabot  and  Lee  1991,  Lilly  1992,  Zang,  Street  and  Koseff  1993,  Vreman,  Geurts  and  Kuerten  1994).  In  this  approach 
constants  which  appear  in  the  base  sgs  models  are  computed  dynamically  as  functions  of  space  and  time  providing  the 
proper  local  amount  of  sgs  scalar  mixing  and  dissipation.  Another  important  idea  in  sgs  modeling  involves  the  use  of 
scale-similarity  laws  which  assume  that  the  largest  of  the  unresolved  scales  have  similar  structure  to  the  smallest  of 
the  resolved  scales.  Among  these  models,  only  that  proposed  by  Bardina,  Ferziger  and  Reynolds  (1983)  satisfies  the 
important  physical  constraint  of  Galilean  invariance  (Speziale  1985). 

Dynamic  modeling  coupled  with  scale  similarity  ideas  have  been  very  useful  in  sgs  modeling  of  nonreacting 
flows.  This  concept  is  currently  being  exploited  in  LES  of  reacting  flows.  Erlebacher,  Hussaini,  Speziale  and  Zang 
(1992)  were  one  of  the  first  to  propose  a compressible  generalization  of  the  dynamic  Smagorinsky  model  and  couple 
this  to  the  scale-similarity  model  to  obtain  the  sgs  mass  and  energy  fluxes.  These  ideas  have  now  been  extended  to 
multicomponent  mixtures  (Oefelein  2001).  With  this  framework  in  place,  the  role  of  combustion  models  are  to  account 
for  the  effect  of  sgs  fluctuations  in  the  thermochemical  variables  and  filtered  chemical  source  terms  over  a wide  range 
of  pressures.  Achieving  this  closure  hinges  on  the  specification  of  accurately  reduced  kinetics  mechanisms  coupled  to 
an  accurate  representation  of  sgs  scalar  mixing  processes. 
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Numerical  Constraints 


It  is  well  known  that  numerical  dissipation  and  dispersion  errors  can  have  significantly  devastating  effects  on  sgs 
models.  The  presence  of  these  errors  depletes  energetic  turbulence  scales  at  the  mid-  to  high-wavenumbers  and  con- 
sequently competes  with  the  models.  When  this  occurs,  the  sgs  models  themselves  often  have  no  effect  on  the  flow, 
and  the  contamination  often  to  leads  erroneous  conclusions.  To  avoid  this  situation,  numerical  methods  used  for  LES 
must  provide  spatially  “non-dissipative”  spectrally  clean  damping  characteristics  out  to  the  smallest  wavenumbers 
coupled  with  simultaneous  local  conservation  of  mass,  momentum  and  total-energy.  Conservation  of  kinetic  energy  is 
particularly  important  when  dynamic  modeling  is  used. 

Co-located  schemes  with  explicit  artificial  dissipation  terms  added  for  stabilization  purposes  have  historically 
failed  to  provide  the  appropriate  spectral  characteristics.  This  is  easily  shown  if  one  compares  the  magnitude  of  the 
residual  associated  with  the  artificial  dissipation  terms  of  a given  scheme  to  that  associated  with  a given  sgs  model.  The 
former  is  always  orders  of  magnitude  greater,  even  for  higher-order  schemes.  Unfortunately,  this  fact  precludes  a wide 
class  of  flow  solvers,  including  the  trivial  conversion  of  most  RANS  based  codes.  Staggered  grid  algorithms  fashioned 
after  the  pioneering  work  of  Harlow  and  Welch  ( 1965),  on  the  other  hand,  have  been  shown  to  give  acceptable  spectral 
characteristics.  Specialized  schemes  of  this  type  are  currently  the  workhorse  of  the  LES  community. 

Because  of  the  obvious  advantages,  grid  stretching  functions  are  widely  used  for  LES  just  as  they  are  for  RANS 
calculations.  One  additional  constraint  associated  with  LES,  however,  is  that  the  energetic  scales  must  be  resolved  on 
grids  that  minimize  commutation  errors.  This  requirement  imposes  strict  grid  stretching  constraints  which  precludes 
the  use  of  typical  RANS  grids.  Instead,  grids  must  be  constructed  with  much  more  restricted  grid  stretching  and  grid 
aspect  ratios.  The  stretching  ratio  associated  with  adjacent  cells  in  a given  coordinate  direction  should  never  exceed 
10  percent,  and  grid  aspect  ratios  greater  than  100  are  rarely  acceptable.  The  issues  outlined  above  represent  mini- 
mal requirements  and  non-adherence  can  lead  to  diminished  broadband  resolution  and  significant  high  wavenumber 
contamination. 


Phenomenological  Case  Studies 

Development  efforts  conducted  by  Oefelein  (2001)  over  the  past  several  years  have  led  to  a massively  parallel  software 
package  which  incorporates  the  general  requirements  for  LES  described  in  the  previous  section.  The  effort  was  driven 
by  two  mutually  dependent  objectives.  The  first  was  to  develop  improved  models  suitable  for  performing  high-fidelity 
LES  of  the  complex  phenomena  described  above.  The  second  was  to  develop  a high-performance  parallel  algorithm 
which  supported  the  implementation  of  large-scale  simulations.  Emphasis  was  placed  on  a general  treatment  of  phe- 
nomenologically complex  reacting  multiphase  flows,  including  the  seven  fundamentally  important  flow  characteristics 
described  above.  These  flow  characteristics  are: 

1.  Dense  near-critical  and  supercritical  fluid  mixture  properties. 

2.  Broadband  turbulent  mixing  over  a wide  range  of  scales. 

3.  High  pressure  chemical  kinetics. 

4.  Strong  multicomponent  property  gradients. 

5.  Dominant  preferential  diffusion  processes. 

6.  Anomalous  multiphase  interfaces. 

7.  Geometrically  dominated  (wall-bounded)  three-dimensional  evolution. 

The  baseline  Eulerian-Lagrangian  framework  solves  the  filtered  conservation  equations  of  mass,  momentum,  total- 
energy  and  species  using  a staggered  grid  methodology  analogous  to  that  pioneered  by  Harlow  and  Welch  (1965) 
in  generalized  curvilinear  coordinates.  Dual-time  stepping  is  used  with  a unified  all  Mach  number  preconditioning 
technique.  The  algorithm  accommodates  fully  implicit  time  advancement  using  a fully  explicit  multistage  scheme  in 
pseudo-time.  This  scheme  exhibits  excellent  parallel  efficiency  and  scalability  attributes.  The  implicit  formulation  is 
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A-stable  which  allows  one  to  set  the  time  step  based  solely  on  accuracy  considerations.  It  accommodates  arbitrary 
equations  of  state,  thermochemical,  thermodynamic  and  transport  processes,  and  provides  full  thermophysical  cou- 
pling over  a wide  range  of  conditions.  It  accommodates  intermediate  complex  geometric  features  while  maintaining 
the  accuracy  of  structured  spatial  stencils.  The  parallel  paradigm  employs  distributed-memory  message-passing  using 
MPI,  the  single-Program-Multiple-Data  (SPMD)  model  and  structured  multiblock  domain  decomposition. 

Following  are  three  examples  which  give  insight  into  the  intricate  nature  of  the  problem  and  the  current  state 
of  the  art  with  respect  to  LES.  The  first  is  an  early  set  of  results  which  illustrate  the  prevalence  of  items  1-7  listed 
above  when  modeling  high  pressure  mixing  and  combustion  in  liquid-oxygen  (LOX),  hydrogen  systems.  These  results 
represent  a first  preliminary  attempt  at  simulating  these  phenomena  and  demonstrate  a capability  to  handle  the  extreme 
complex  thermophysical  flow  characteristics.  The  second  example  is  an  LES  of  a low-Mach-number,  high-Reynolds- 
number,  particle-laden  channel  flow.  This  is  an  extremely  difficult  case  to  handle  with  a compressible  flow  solver.  The 
results  demonstrate  a capability  to  handle  these  extremes,  the  ability  of  LES  to  capture  transient  broadband  turbulent 
mixing  over  a wide  range  of  scales,  and  also  the  ability  to  resolve  the  geometrically  dominated  (wall -bounded)  three- 
dimensional  evolution  of  the  flow.  The  last  set  of  results  also  demonstrates  these  advantages  in  a more  complex 
geometry,  and  additionally  the  ability  of  LES  to  simulate  spray  characteristics  and  particle  dispersion. 

High-Pressure  Mixing  and  Combustion  in  LOX-H2  Systems 

Figure  8 shows  contours  of  density,  temperature,  and  FU,  O2,  OH  and  H9O  mass  fractions  in  the  near-field  injector 
region  of  a hydrogen-liquid-oxygen  shear  layer.  The  two  streams  are  separated  by  a 0.5  mm  LOX  post.  The  pres- 
sure is  10.1  MPa  (100  atm).  The  hydrogen  (upper  stream)  and  oxygen  (lower  stream)  velocities  are  125  and  30  m/s, 
respectively.  The  injection  temperatures  are  150  K and  100  K . These  early  calculations  were  performed  using  the 
theoretical-numerical  framework  developed  by  Oefelein  (1997)  and  represent  a first  attempt  at  simulating  such  pro- 
cesses. The  matrix  of  conditions  considered  were  fashioned  after  the  flow  visualization  studies  conducted  by  Mayer 
and  Tamura  (1996).  Emphasis  was  placed  on  the  near-field  flow  processes  in  the  vicinity  of  the  post.  The  conditions 
selected  produce  a supercritical  hydrogen  stream  and  a liquid-oxygen  stream  which  undergoes  a transcritical  change 
of  state  within  the  mixing  layer.  Inlet  velocity  profiles  were  generated  assuming  fully  developed  turbulent  flow  and 
a heat  conduction  model  was  applied  to  the  splitter  plate  to  provide  a realistic  energy  flux  distribution  at  the  walls. 
Nonreflecting  outflow  conditions  were  imposed  at  the  exit  and  inviscid,  adiabatic,  and  noncatalytic  conditions  were 
imposed  at  the  transverse  boundaries. 

These  results  illustrate  the  prevalence  of  items  1-7  listed  above.  Transcritical  mixing  induces  a vortical  structure 
within  the  injected  hydrogen  stream  which  is  analogous  to  that  produced  by  a backward  facing  step.  This  structure 
emanates  from  the  boundary  layer  upstream  of  the  post  and  is  amplified  by  interactions  within  the  shear  layer  and 
coalescence  downstream  with  adjacent  vortices.  The  oxygen  stream,  on  the  other  hand,  proceeds  unimpeded  in  an 
essentially  straight  line.  Because  of  the  liquid-like  characteristics  of  the  oxygen  stream,  an  extremely  large  density 
gradient  exists  within  this  region.  Note  that  the  change  in  density  is  on  the  order  of  1000  to  1.  Diminished  mass 
diffusion  rates  are  also  evident.  The  combined  effect  produces  a fuel  rich  flame  which  anchors  itself  to  the  oxygen  jet 
and  behaves  in  a qualitatively  similar  manner  as  the  diffusion  dominated  flame  depicted  in  Figs.  2 and  3.  Combustion 
occurs  at  near  stoichiometric  conditions  and  produces  a wake  which  effectively  separates  the  hydrogen  and  oxygen 
streams  as  the  flow  evolves  downstream.  Results  highlight  the  effect  of  the  momentum  flux  ratio  on  flame-holding 
dynamics,  the  dominating  effect  of  the  density  gradient,  and  the  impact  of  diminished  mass  diffusion  rates  which 
accompany  the  liquid-like  behavior  of  near-critical  fluids. 

Low-Mach-Number,  High-Reynolds-Number  (Particle-Laden)  Channel  Flow 

Three  particularly  relevant  effects  induced  by  interphase  coupling  are  turbulence  modulation  which  involves  the  damp- 
ing of  gas  phase  turbulence  by  particulates  accommodating  to  turbulent  motion,  turbulence  generation  which  involves 
the  production  of  gas  phase  turbulence  due  to  the  presence  of  particle  wakes,  and  liquid  deformation  and  breakup 
processes  and  the  resultant  effect  on  interphase  exchange  processes.  As  part  of  an  effort  to  treat  these  phenomena 
systematically  a series  of  LES  calculations  are  being  performed  using  the  algorithmic  framework  describe  above  and 
compared  to  the  experimental  data  acquired  by  Kulick,  Fessler  and  Eaton  (1994).  These  experiments  characterize  the 
interactions  between  various  particle  loading  conditions  and  the  fluid  turbulence  in  the  well-defined  confines  of  a tur- 
bulent channel.  Particles  are  selected  to  respond  to  some,  but  not  all  scales  of  turbulent  motions.  Gas  phase  velocities 
were  measured  to  investigate  the  means  by  which  particles  attenuate  turbulence. 
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Figure  8:  Contours  of  density,  temperature,  and  H2,  O2,  OH  and  HoO  mass  fractions  in  the  near-field  injector  region. 
Chamber  pressure  is  10.1  MPa  ( 1 00  atm),  hydrogen  (upper  stream)  and  oxygen  (lower  stream)  velocities  are  1 25  and 
30  m/s,  respectively,  and  injection  temperatures  are  150  K and  100  K. 
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6 3 


Figure  9:  Instantaneous  contours  of  the  streamwise  component  of  vorticity  in  a high  Reynolds  number,  particle  laden 
channel.  Half-height  is  h = 20  mm , mean  centerline  velocity  is  Uc l = 10.5  m/s , Reynolds  number  based  on  h is 
Re^  = 13800,  friction  velocity  is  uT  = 0.49  m/s  and  Reynolds  number  based  on  uT  is  Rer  = 645. 
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Figure  10:  Velocity  vectors  and  contours  of  the  streamwise  component  of  vorticity  in  the  y ■ 
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Figure  11:  Time-averaged  streamwise  component  of  velocity  compared  with  classic  law-of-the-wall  profiles. 


Figure  9 shows  the  instantaneous  contours  of  the  streamwise  component  of  vorticity  in  the  turbulent  low-Mach- 
number,  high-Reynolds-number,  particle  laden  channel  of  Kulick  et  al.  (1994).  The  channel  half-height  is  h = 20  mm. 
The  mean  centerline  velocity  is  Uc\  = 10.5  m/s  and  the  Reynolds  number  based  on  h is  Re  h = 13800.  The  friction 
velocity  is  uT  = 0.49  m/s  and  Reynolds  number  based  on  ur  is  Rer  = 645.  The  domain  dimensions  in  the  stream- 
wise,  transverse  (wall-normal),  and  spanwise  directions,  respectively  are  Qh  x 2h  x 3 h.  The  primary  grid  is  composed 
of  1003  hexahedral  cells. 

This  second  set  of  results  demonstrate  the  ability  of  LES  to  capture  transient  broadband  turbulent  mixing  over 
a wide  range  of  scales  and  the  ability  to  resolve  the  geometrically  dominated  three-dimensional  evolution  of  turbu- 
lent flows.  It  also  demonstrates  the  ability  of  the  algorithmic  framework  to  handle  low-Mach-number  flows  in  the 
incompressible  limit.  The  detailed  broadband  structure  associated  with  the  baseline  conditions  described  above  are 
shown  in  Fig.  10.  Here  velocity  vectors  are  plotted  with  contours  of  the  streamwise  component  of  vorticity  in  the 
y — z plane.  Figures  1 1 and  12  show  the  time-averaged  streamwise  component  of  velocity  compared  with  classic 
law-of-the-wall  profiles  and  the  streamwise  and  wall-normal  root-mean-square  components  of  velocity  compared  with 
experimental  data  points  obtained  by  Kulick  et  al.  (1994).  These  results  are  typical  when  dynamic  modeling  is  used 
with  an  appropriately  sized  grid. 

Swirling  Particle-Laden  Flow  in  Coannular  Dump-Combustor 

In  addition  to  broadband  turbulence  structure,  this  last  set  of  results  demonstrates  the  ability  of  LES  to  simulate  spray 
characteristics  and  particle  dispersion.  Obtaining  high-fidelity  solutions  of  reacting  sprays  hinges  on  the  application 
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Figure  12:  Streamwise  and  wall-normal  root-mean-square  components  of  velocity  compared  with  experimental  data 
points  obtained  by  Kulick  et  al.  (1994). 
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Figure  13:  Cross-section  of  the  computational  domain  and  grid  employed. 


of  methods  and  models  which  accurately  describe  momentum  coupling  and  subgrid-scale  modulation  of  turbulence, 
mass  and  energy  coupling  and  subgrid-scale  scalar  mixing,  and  the  combustion  processes  induced  as  a consequence. 
As  part  of  an  effort  to  treat  these  phenomena  systematically  a series  of  LES  calculations  have  been  performed  and 
compared  to  the  experimental  data  acquired  by  Sommerfeld  et  al.  (Sommerfeld  and  Qiu  1991,  Sommerfeld,  Ando 
and  Wennerberg  1992,  Sommerfeld  and  Qiu  1993).  These  experiments  characterize  a swirling  particle-laden  flow  in  a 
model  coannular  combustion  chamber  and  effectively  isolate  the  effects  related  to  momentum  coupling.  The  primary 
objectives  here  were  to  gain  a clearer  understanding  of  the  effectiveness  and  feasibility  of  current  models  and  to  gain 
a quantitative  understanding  of  potential  model  limitations  by  analyzing  the  characteristic  fluid  dynamic  scales  of 
importance. 

Sommerfeld  et  al.  provides  detailed  measurements  of  swirling  particle-laden  flow  in  a chamber  which  consists 
of  a sudden  pipe  expansion  with  a centered  (primary)  and  annular  (secondary)  jet  discharging  into  a cylindrical  test 
section.  A cross-section  of  the  computational  domain  is  given  in  Fig.  13.  The  region  of  interest  is  shown  in  Fig.  14. 
The  primary  jet  has  a radius  of  r/R  — 0.5  and  is  laden  with  glass  beads  with  a mean  particle  diameter  of  45  jim 
distributed  between  20  and  80  fim.  The  secondary  jet  extends  over  a radial  interval  of  0.59  < r/R  < 1 and  is  injected 
with  a swirling  azimuthal  velocity  component.  The  relevant  flow  conditions  and  particle  properties  are  summarized 
in  Table  1.  Particles  are  injected  in  the  primary  jet  in  equilibrium  with  the  gas-phase  in  a manner  that  matches  the 
experimental  distribution.  A series  of  one-component  phase-Doppler  anemometer  measurements  were  made  along 
cross-sections  at  the  8 axial  locations  indicated.  Gas-phase  and  particle -phase  mean  and  rms  velocity  components 
were  acquired  along  with  simultaneous  measurements  of  the  particle  size  and  mass  flux  distributions. 

Figure  15  shows  1 mm  thick  cross-sections  of  the  instantaneous  particle  distributions  for  Case  2 superimposed 
on  the  corresponding  turbulent  velocity  field.  For  plotting  purposes,  the  superimposed  particle  distribution  shown  is 
quite  thick  relative  to  the  particle  size  distribution.  One  should  not  infer  from  this  figure  that  collisions  are  important. 
Analysis  of  the  particle  number  densities  throughout  the  flowfield  indicate  that  they  are  not.  At  any  instant  in  time  there 
are  approximately  2.5  million  particles  being  tracked,  with  two-way  coupling  applied  between  the  gas  and  particles. 
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Figure  14:  Region  of  interest  in  the  Sommerfeld  configuration. 


Table  1:  Flow  conditions  and  particle  properties  used  in  the  Sommerfeld  experiments. 


Case  1 

Case  2 

Gas  Phase  (Air): 

Flow  rate  in  primary  jet  g/s 

9.9 

6.0 

Flow  rate  in  secondary  jet,  g/s 

38.3 

44.6 

Inlet  Reynolds  numbera 

26200 

27250 

Swirl  number 

0.47 

0.49 

Temperature,  K 

300 

Particle  Phase: 

Loading  ratio  in  primary  jet 

0.034 

0. 17  6 

Flow  rate,  g/s 

0.34 

1.0 

Mean  diameter,  pm 

45.5 

Density  ratio,  pP/pj 

2152 

a Based  on  total  volume  flow  rate. 
6 5 x Case  1. 
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Figure  15:  One  millimeter  cross-sections  of  the  instantaneous  particle  distribution  for  Case  2 superimposed  on  the 
corresponding  turbulent  velocity  field  ( Rref  — 32  mm,  C7re/  = 12.9  m/s). 


Tracking  this  number  of  particles  is  significant  since  it  verifies  the  feasibility  of  employing  large  numbers  of  physical 
particles  and  eliminates  the  need  to  implement  classical  “parcel”  approximations. 

The  mean  flow  characteristics  for  Cases  1 and  2 are  shown  in  Fig.  16.  This  figure  shows  the  time-averaged, 
azimuthally-averaged  gas-phase  velocity  field.  Key  features  of  the  flow  include  primary  and  secondary  recirculation 
zones,  a stagnation  point  in  the  core  region  and  a reattachment  point  on  the  outer  wall.  The  location  of  these  points 
coincide  with  measured  results  to  within  5 % for  both  cases.  Figure  17  shows  representative  comparisons  of  (a)  the 
mean  axial  gas  velocity,  (b)  the  corresponding  mean  axial  particle  velocity,  (c)  the  mean  particle  diameter,  and  (d) 
the  mean  particle  momentum  flux  for  Case  2 at  respective  axial  stations.  The  agreement  between  the  measured  and 
calculated  results  is  excellent  and  similar  agreements  have  been  obtained  with  respect  to  the  entire  experimental  data 
set. 

After  validating  the  LES  methodology  with  the  Sommerfeld  data  subsequent  calculations  were  performed  within 
the  same  configuration  to  gain  a more  quantitative  understanding  of  the  relevant  modeling  issues.  Here,  the  validated 
broadband  characteristics  inherent  to  the  LES  methodology  were  used  to  obtain  additional  data  that  could  not  be 


Figure  16:  Mean  flow  gas-phase  stream  function  and  velocity  vectors. 
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Figure  17:  Comparisons  of  measured  (o)  and  calculated  (-)  time-averaged  profiles  of  (a)  the  mean  axial  gas  velocity, 
(b)  the  corresponding  mean  axial  particle  velocity,  (c)  the  mean  particle  diameter,  and  (d)  the  mean  particle  momentum 
flux  for  Case  2. 
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Figure  18:  Mean  particle  Reynolds  number  corresponding  to  Case  2:  (—)  all  classes;  ( ■)  dp  = 30  ± 5 pm;  (— ) 
dp  = 45  ± 5 pm;  (— ) dp  = 60  ± b pm . 


acquired  from  the  experiment.  In  addition  to  the  enhanced  quantitative  insight,  this  data  was  used  to  determine  the 
envelope  of  conditions  to  be  characterized  by  the  model  development  effort.  Figure  18  presents  a representative  set 
of  results.  Here  profiles  of  the  mean  particle  Reynolds  number  for  Case  2 are  given  at  respective  axial  stations  for  all 
size  classes,  d = 30  ± 5 pm,  d = 45  ± 5 pm,  and  d = 60  ± 5 pm.  Other  key  scales  and  parameters  obtained  include 
1)  particle  characteristics  such  as  the  particle  number  density,  volume  fraction,  Sauter  mean  diameter  etc.,  2)  flow 
characteristics  such  as  the  turbulence  intensity,  Kolmogorov  microscales  and  Taylor  microscales,  and  3)  time-scales 
such  as  the  Stokes  number  and  particle  relaxation  time. 


Conclusions 

Improvements  in  computer  speed  and  capacity  have  made  application  of  the  large  eddy  simulation  technique  feasi- 
ble for  increasingly  complex  flows.  The  discussion  and  results  presented  have  outlined  the  key  phenomenological 
trends  and  have  demonstrated  model  performance  and  accuracy  requirements.  Results  have  also  highlighted  various 
intricacies  associated  with  transcritical  and  supercritical  phenomena,  highlighted  the  effect  of  pressure  on  near-critical 
mixing  and  combustion  processes,  and  provided  increased  insights  into  the  theoretical  and  numerical  methodolo- 
gies employed.  Currently,  coupled  turbulent  mixing  and  the  treatment  of  thermodynamic  and  transport  processes  in 
laboratory-scale  geometries  can  be  handled  quite  accurately.  The  treatment  of  dilute  spray  dispersion  and  vaporization 
processes  are  also  handled  well.  Accurate  treatment  of  turbulent  premixed  and  non-premixed  flame  phenomena,  on  the 
other  hand,  is  still  pending,  as  is  the  treatment  of  atomization  processes  and  interface  dynamics.  The  latter  still  lacks 
a strong  theoretical  basis.  Current  efforts  are  focused  on  model  assessment  and  validation  at  realistic  device-scale 
conditions  and  analysis  of  validated  systems  to  systematically  characterize  the  relevant  time-scales,  length-scales,  and 
other  key  parameters  which  are  of  direct  importance  to  sgs  model  development. 
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